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Abstract 



We derive the equations of motion for binary systems of compact bod- 
ies in the post-Newtonian (PN) approximation to general relativity. Results 
are given through 2PN order (order (f/c) 4 beyond Newtonian theory), and 
for gravitational radiation reaction effects at 2.5PN and 3.5PN orders. The 
method is based on a framework for direct integration of the relaxed Ein- 
stein equations (DIRE) developed earlier, in which the equations of motion 
through 3.5PN order can be expressed in terms of Poisson-like potentials that 
are generalizations of the instantaneous Newtonian gravitational potential, 
and in terms of multipole moments of the system and their time derivatives. 
All potentials are well defined and free of divergences associated with inte- 
grating quantities over all space. Using a model of the bodies as spherical, 
non-rotating fluid balls whose characteristic size s is small compared to the 
bodies' separation r, we develop a method for carefully extracting only terms 
that are independent of the parameter s, thereby ignoring tidal interactions, 
spin effects, and internal self-gravity effects. Through 2.5PN order, the re- 
sulting equations agree completely with those obtained by other methods; the 
new 3.5PN back-reaction results are shown to be consistent with the loss of 
energy and angular momentum via radiation to infinity. 

04.30.-w, 04.25.Nx 



1 



Typeset using REVTpjX 



I. INTRODUCTION AND SUMMARY 



This is the second in a series of papers which will treat motion and gravitational radiation 
in the post-Newtonian approximation to general relativity. While this is a problem that dates 
back to the beginnings of general relativity, it has recently taken on added observational 
importance because of the need for extremely accurate theoretical gravitational waveform 
templates for analysis of data taken by laser interferometric gravitational-wave detectors [|T| . 
Specifically, for waves from inspiralling binary systems of compact objects (neutron stars or 
black holes), equations of motion and gravitational waveforms accurate to at least third post- 
Newtonian order (order (v/c) 6 ) beyond the initial Newtonian or quadrupole approximation 
are needed. 

In paper I 0, we laid out the foundations of our method of Direct Integration of the 
Relaxed Einstein Equations (DIRE). We rewrote the Einstein equations as a flat spacetime 
wave equation together with a harmonic gauge condition (the "relaxed" Einstein equations), 
and solved them formally in terms of a retarded integral over the past null cone of the 
field point. Because the "source" contains both the material stress-energy tensor and the 
stress-energy contributions of the gravitational fields themselves, it was necessary to iterate 
the integrals repeatedly to obtain successively higher-order approximations to a solution in 
powers of e ~ (v/c) 2 ~ (Gm/rc 2 ). Each power of e represents one "post-Newtonian" (PN) 
order in the series (e 1 / 2 represents one half, or 0.5PN orders). Despite the fact that the field 
contributions to the integrals extend over all spacetime, we demonstrated that no infinite 
or ill-defined integrals occurred, even in slow-motion, multipole expansions, and found a 
simple prescription for evaluating the finite contributions of all integrals. This was true for 
calculations of the metric both in the near zone and in the far zone. 

To complete the solution of Einstein's equations, one needs equations of motion for the 
system. For this, one needs the spacetime metric evaluated for field points within the near 
zone, corresponding to a sphere of radius 1Z ~ one gravitational wavelength. In Paper 
I, we expressed this near-zone metric explicitly through order e 7//2 beyond the Newtonian 
approximation, corresponding to 3.5 post-Newtonian (PN) order, in terms of instantaneous, 
Poisson-like integrals and their generalizations, of the form, for example, 

P<J) = -j-f f^d 3 x', (1.1) 
Arc Jm |x — x'| 

where the integration is confined to the near zone Ai, and only that part of the integral 
that is independent of 1Z is kept. 

It is the purpose of this paper to evaluate these integrals explicitly for a binary system of 
non-spinning, spherically symmetric bodies whose size is much smaller than their separation. 
We will carry this evaluation through 2PN order, and will also evaluate the leading radiation- 
reaction contributions at 2.5PN order, together with the first post-Newtonian corrections to 
radiation reaction, at 3.5PN order. The extremely lengthy derivation of the non-radiative 
3PN contributions will be reserved for future publications. The resulting equations have the 
form 

g? 2 x 771 in 

-jjj = + ^[ n (Ap N + Mpn + A 3PN ) + fv(B PN + B 2PN + B 3PN )} 

8 Tfl Tfl 

+TV— — [ rn ( A 2.5PN + A 3 . 5PN ) - v(B 2 . 5PN + B 3 , 5PN )] , (1.2) 

T T 
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where x = Xi — X2, r = |x|, n = x/r, m = mi + rn.2, 77 = mim2/m 2 , v = vi — V2, and 
r = dr/dt. We use units in which G = c = 1. The leading term is Newtonian gravity. The 
other terms on the first line are the "conservative" or non-dissipative terms, of even PN 
order, while those on the second line are dissipative radiation-reaction terms, of odd-half 
PN order. The coefficients A and B are given explicitly by 

3 

A PN = -(1 + 3t])v 2 + -7]f 2 + 2(2 + r])m/r , 

B PN = 2(2 - 77) , (1.3a) 

1 3 
A2PN = —^(3 — 4r/)t> 4 + 2^(13 ~~ 4:T])v 2 m/r + —77(3 — Ar])v 2 r 2 

15 3 

+ (2 + 25?/ + 2rj 2 )r 2 m/r 77(1 - 3r?)r 4 - -(12 + 29r/)(m/r) 2 , 

8 4 

1 3 1 

5 2 pv = 2^ 15 + 4r ^ 2 ~ 2^ 3 + 2r/ ^ 2 ~ 2 ^ + 4l7/ + 8r/2 ) m / r ' ^ L3b ) 

^2.5pjv = 3w 2 + —m/r, 

B2.5PN = v 2 + 3m/r , (1.3c) 

^3.5PV = -4( 61 + 70r /) w4 - i( 519 - 1267r/)w 2 m/r + —(19 + 2r])v 2 r 2 
28 42 4 

1 23 
--(147 + 188r])r 2 m/r - 70f 4 - —(43 + Ur])(m/r) 2 , 

5 3 . 5 pv = "4( 313 + 42r ^)^ 4 + i( 205 + 777r])v 2 m/r + ^(113 + 2r/)w 2 r 2 
28 42 4 

- — (205 + 424r/)r 2 m/r - 75r 4 - —(1325 + 546r?)(m/r) 2 . (1.3d) 
1 2 42 



The 1PN coefficients are standard; the 2PN and 2.5PN terms agree completely with 
results derived by Damour and Deruelle [||,|J , Kopeiken and Grishchuk |||6| , Blanchet et al. 
0, and Itoh et al. ||. However, the 2.5PN terms differ from those derived, for example from 
a Burke-Thorne type radiation-reaction potential given by §rr = —(l/5)x' l x : 'd 5 I <t: ' > /dt 5 , 
where X <4:,> is the system's trace-free quadrupole moment (see, e.g. §36.11 of ||). Iyer and 
Will JTTJJTT]] showed that there is a two-parameter "gauge" freedom in the radiation reaction 
equations at 2.5PN order (and a six-parameter freedom at 3.5 PN order), within which the 
different equations of motion yield identical results for the net energy and angular momentum 
radiated. The 2.5PN terms shown are thus observationally equivalent to the Burke-Thorne 
radiation reaction equations. The 3.5PN terms are new; it can be shown that they correspond 
to a specific choice of the six Iyer- Will gauge parameters, and thus automatically generate 
the proper post-Newtonian corrections to energy and angular momentum loss. Schafer and 



Jaranowski [IT2HT3] also derived 2PN, 2.5PN and 3.5PN contributions to the equations of 



motion in a Hamiltonian formulation. The 3PN contributions to the equations of motion 
have also been reported by several groups |fT5|-r9"|. 

The remainder of this paper is devoted to the details supporting these results. In Sec. 
[TI] we review the basic equations needed to find equations of motion to the order needed. 
Section [111] specializes to binary systems of spherical "pointlike" bodies, and derives the 
equations of motion for each body correct to 2PN order. Section [IV| repeats the process 
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for the 2.5PN and 3.5PN terms. In Sec. [V| we transform to an effective one-body relative 



equation of motion. Concluding remarks are made in Sec. [VI[ Detailed points are reserved 
for a series of Appendices. 

Our conventions and notation generally follow those of ||,^0|. Greek indices run over 



four spacetime values 0, 1, 2, 3, while Latin indices run over three spatial values 1, 2, 
3; commas denote partial derivatives with respect to a chosen coordinate system, while 
semicolons denote covariant derivatives; repeated indices are summed over; rf = rj^ = 
diag(-l, 1, 1, 1); g = det(^); = (a« + a^)/2; = (a ij - a ji )/2; e ijk is the totally 
antisymmetric Levi-Civita symbol (e 123 = +1). We use a multi-index notation for products 
of vector components and partial derivatives, and for multiple spatial indices: x*i— k = 
x l x 3 ' . . . x k , V lJ '" fc = V l V J . ..V fc , with a capital letter superscript denoting an abstract 
product of that dimensionality: x® = x ll x 12 ...x lq and = V n V l2 ...W Also, for a tensor 
of rank Q, f® = f lll ^- l i . The notation f®g® denotes a complete contraction over the q 
indices. Spatial indices are freely raised and lowered with and Sij. 

II. EQUATIONS OF MOTION OF COMPACT BINARY SYSTEMS: BASIC 

EQUATIONS 

A. Structure of the near-zone metric to 3.5PN order 

We begin by reviewing key results from Paper I. We defined the "field" h af3 by 

h ^ = v aP_ { _ g) l/2 g a^ (21) 

where g a @ is the space time metric. In deDonder or harmonic coordinates defined by the 
gauge condition h a ^ \p = 0, the Einstein equations take the form 

Uh al3 = -167TT af3 , (2.2) 

where □ is the flat-spacetime wave operator, and r Q/3 is made up of the material stress- 
energy tensor T 0/3 and the contribution of all the non-linear terms in Einstein's equations. 
We defined a notation for specific components of the field h af3 : 

N = h 00 ~ 0(e) , 
K l = h 0i ~ 0(e 3 / 2 ) , 
B ij = h ij ~ 0(e 2 ) , 

B = Y,h u ^ 0(e 2 ) , (2.3) 

i 

where we show the leading-order dependence on e in the near zone. To the necessary orders 
for calculating the equations of motion to 3.5PN order, the components of the physical 
metric are given in terms of N, K l , B^ and B by 

#oo = -(1 ~-N+ -N 2 - -N 3 + — iV 4 ) + -B(l --N+ -N 2 ) 
yuu v 2 8 16 128 ' 2 v 2 8 ; 

+ im»B ij - -B 2 ) + -K j K j - -NK j K j + 0(e 5 ) , 
4 V 2 ; 2 4 y J 
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g oi = -K\l -In~~B + *-N 2 ) - K*B» + 0(e 9 / 2 ) , 

9H = + \n- l -N 2 + ±N 3 - X -NB + l -K k K k ) 

+B ij - ^B5 lj - K l K j + ]^NB ij + 0(e 4 ) , 
(- g ) = i + N - B - NB + K l K { + 0(e 4 ) . (2.4) 
The potentials N, K\ B^ and B must also be expanded to an appropriate order: 
N = e(4U + eiVx + e 3/2 AT L5 + e 2 N 2 + e 5/2 N 2 . 5 + e 3 N 3 + e 7/2 N 3 . 5 ) + 0(e 5 ) , 

K i = e 3/2( 4V H + ^ + e 3/2^ 5 + ^ + e V2 X i ^ + ( £ 9/2) ? 
B = e 2 (jBl + e l/2 jBl 5 + eB2 + ^2^ 5 + + e 5/ 2j g 3 s) + 0(e 5 } ^ 

B ij = e 2 (£^ + e l/2 B% + e^' + e 3/2 ^ 5 ) + 0(e 4 ) , (2.5) 

where the subscript on each term indicates the level (1PN, 2PN, 2.5PN, etc.) of its leading 
contribution to the equations of motion. Notice that our separate treatment of B and B 1 ^ 
leads to the slightly awkward notational circumstance that, for example, B™ = B\. In Paper 
I, we obtained explicit near-zone expressions for each of the terms in Eq. (|2.5| ) in terms of 
Poisson-like potentials 



p(f)=^f f^4 rfV ' v 2 p(/)=-/ 

47T Jm x — x' 



I a(t x-)/(t x0 ^ / = 

JM |X — X'| 

and various generalizations [for definitions, see Paper I, Eqs. (4.10) - (4.16)], and in terms 
of source multipole moments, such as X y = J M r 00 x^d 3 x [for definitions, see Paper I, Eqs. 
(2.14) and (4.5) - (4.7)]. The integrations are over a constant time hypersurface M. that 
extends to a radius 1Z ~ one gravitational wavelength from the source. In Paper I, we showed 
that, even for Poisson potentials where the function / does not have compact support, 
all contributions to the field h^ u from the integration over M. that depend on 1Z cancel 
corresponding contributions from that part of the field point's past null cone that is outside 
TZ, and thus that any ^-dependent terms that appear in a given integral can simply be 
discarded. 

The potentials given in Paper I were expressed in terms of specific source densities given 

by 

a = T 00 + T u , 
a* = T 0i , 

a ij = T ij . (2.7) 



Ua = f ^2^1d 3 x' = P(4vra) = S(l) , 

JM |x — x'| 

V i = f <^y} d 3 x > = p( 47nT i) = gim . (2.8) 

JM X — X' 



For example, in Eq. (|2.5| 
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Explicit expressions for the remaining terms in Eq. ( |2.5|) can be found in Paper I, Eqs. 
(5.2), (5.4), (5.8), (5.10), (6.2), and (6.4). 



B. Model of the material sources 

We model the material sources in the binary system as perfect fluid, having stress-energy 
tensor 

T afB = (p + p)u a u p + pg aP , (2.9) 

where p and p are the locally measured energy density and pressure, respectively, and u a 
is the four-velocity of an element of fluid. We will assume the bodies to be non-rotating 
(the effects of spin will be treated in future publications), spherically symmetric in their 
comoving rest frames, and small compared to their separation, so that tidal distortions can 
be ignored. 

Our goal is to determine all contributions to the equations of motion that are independent 
of the internal structure, size, and shape of the bodies. We are less interested in formal rigor 
than in having a robust method that captures all the effects without missing any. One 
approach to this has been to assume a "delta-function" or distributional form for the stress- 
energy tensor. This has been criticized because such a source is fundamentally incompatible 
with general relativity, and because it leads to divergences related to the infinite self-field of 
a point mass. A number of methods have been developed in order to extract the finite part 
of such divergent expressions, including the Hadamard partie finie technique (for a recent 
review, see ^TJ. Another approach is related to that of Einstein, Infeld and Hoffman (EIH) 
p2|| : expand the vacuum Einstein equations in a post-Newtonian expansion and match the 
solutions to fields representing the near-zone, Schwarzschild-like field of a static, spherical 
body. The consistency conditions imposed by the matching lead to constraints on the motion 
of the bodies that yield the equations of motion. This has been carried out to 2.5PN order 
by Itoh et al. P^,P1, exploiting a "strong- field point particle" scaling method of Futamase 

A third approach is to treat the bodies realistically as fluid balls with internal energy, 
supported against their self gravity by pressure governed by an equation of state. In this 
case, the mass of each body is composed of rest mass, internal energy and self-gravitational 
binding energy, and the center of mass is defined accordingly. However, at Newtonian and 
1PN order, when finite-size effects such as tidal interactions are ignored, it turns out that all 
vestiges of the internal structure are "effaced" , in the language of Damour |4]] , and the final 
1PN equations of motion depend on one and only one mass as defined above. This procedure 
can be seen in detail, for example, in [f2^], §6.2, where the calculation is actually carried out 
in the parametrized post-Newtonian (PPN) framework, which encompasses a class of metric 
theories of gravity. In many alternative theories, such as scalar-tensor gravity, the effacement 
is violated, and the equations of motion depend on various masses, such as inertial mass 
m, active gravitational mass tua, and passive gravitational mass mp, which may differ by 
amounts depending on the bodies' gravitational binding energy. In GR, the PPN parameters 
are such that all three masses are identical. 

In fact, the 1PN equations of motion derived from this method are identical to those 
obtained from a "delta" function method in which one systematically throws away all terms 
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that are singular when evaluated on each body's world line. At 1PN order, the results are in 
keeping with the idea that general relativity satisfies the Strong Equivalence Principle (see 
§3.3 of |[25||), part of which implies that the motion of bound bodies is independent of their 
internal structure, provided that tidal effects can be ignored. Kopeikin || extended this to 
2PN order, with results consistent with the Strong Equivalence Principle. 

Our approach will be intermediate between the "delta-function" model and the full equi- 
librium fluid ball method. We will neglect pressure p and internal energy density, and treat 
the bodies as balls of baryons characterized by the "conserved" baryon mass density p*, 
given by 

p* = mn^gu°, (2.10) 

where m is the rest mass per baryon, n is the baryon number density, and g = det(g M ;,). 
From the conservation of baryon number, expressed in covariant terms by {nu a ). a = = 
(y/—gnu a ) ya , we see that p* obeys the non-covariant, but exact, continuity equation 

dp*/dt + V • (p*v) = , (2.11) 

where v l = u l /u°, and spatial gradients and dot products use a Cartesian metric. In terms 
of p*, the stress-energy tensor takes the form 

T afi = p\-g)- 1,2 u% a v^ (2.12) 

where v a = (l,v % ). We define the baryon rest mass, center of baryonic mass, velocity and 
acceleration of each body by the formulae 

m A = / p*d 3 x , 

J A 

x A = (l/m A ) ( p*xM 3 £, 

J A 

va = dx A /dt = (l/m A ) / p*v l d 3 x , 

J A 

a A = dv A /dt = (l/m A ) I p*aV 3 x, (2.13) 

J A 

where we have used the general fact, implied by the equation of continuity for p*, that 



0_ 

Of 



J p*(t, x')/(t, x, x')dV = / p*(t, x') ^ + v' • V') /(*, x, x')dV . (2.14) 



C. Structure of the equations of motion to 3.5 PN order 



The definition of the stress-energy tensor in terms of p*, Eq. ( |2.12| ), together with the 
equation of continuity, Eq. (|2.11|) , and the fundamental equations of motion, T a/3 ;/3 = can 
be shown to be equivalent to the geodesic equation u^u^ = for each fluid element. In 
terms of ordinary velocity v l = dx 1 jdt and harmonic coordinate time t, the geodesic equation 
takes the form 
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a 1 = dtf/dt = -I^uV + rJ^uW , (2.15) 

where rig are Christoffel symbols computed from the metric. According to our definitions 
of the baryonic center of mass, velocity and acceleration of each body, we can write the 
coordinate acceleration of the v4-th body in the form 

a\ = (l/m A ) J a p*(-FipV V + Flpv^v^x . (2.16) 

Our task therefore, is to determine the Christoffel symbols through a PN order sufficient for 
equations of motion valid through 3.5PN order using the 3.5PN accurate expressions of the 
metric in Paper I (different components of T^o are need to different accuracy, depending on 
the number of factors of velocity which multiply them); re-express the Poisson potentials 
contained in the metric in terms of p*, rather than in terms of the "densities" a, a 1 and a % \ 
substitute into Eq. (|2.16| ), and integrate over the A-th body, keeping only terms that do not 
depend on the bodies' finite size. 



D. Christoffel Symbols to 3.5PN order 

The fundamental definition 

= \g aP {gpw + - <w) , (2.17) 



together with the form of the metric Eq. (2^4) and the expansions of Eq. ( |2.5|) give the 



Christoffel symbols expanded to the required order: 

Too = ~eU* ~ e 2 + #i) - W a U - W*uA - e^N^ - e 3 Q(iV 2 + B 2 ) - AW CT 

-U a Ni - V:(N( + B() - K\Ut - Wyl + lQU 2 U a + 32U„ViU$) 
-e 7/2 Q(iV 2 . 5 + B 2 . 5 ) - NxtU. - UvNx* - K l 2 . 5 U^ + 0(e 4 ) , (2.18a) 

F° 0l = -dJt - e 2 (±{N? + B() - 4U a U?) - 6 3 (i(J\tf + Bi) - N x Ui - U a N} 

+ W'V? + IQU 2 U^ - e 7 / 2 (±(Nl 5 + BQ - NutU?) + 0(e 4 ) , (2.18b) 

r° = e 2 {W^ + (fjv) + e 3 (ifa - B x )5 lj + ^B% + Kf j) - 16V^ + 4V a k U> k 5 ij ^) 

+t 7/2 Qb& + - ijVi.5**) + 0(e 4 ) , (2.18c) 

r^ = -U£-e (±(N( + Bt) + W; - 8U a uty - e 2 + Bi) +K*- 
-2U C N? - U a B( - m/l ■■;! - W$J a - l6U a V; + SViVi* + 4&1%U* 
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-t?(~{Ng + B£) +Kl- 2N 2 U£ - U a (2Ng + B$) - AU a K l 2 - K\tl a 

- \n x B( - V l a {Ni + B x ) - AN X V: + 2K 2 i VJ' i + 2V'K$ i + 4V'B$ 

+~B 3 2 k B 3 2 k > i - B$U* - ^B${B? + N[ j ) + 12[/ 2 jVf + AU 2 a B( + 2AU (T N l U£ 

+8U a UJB l 2 J + MU*Vi + 32U a U a V: - QAU^V^ - Z2V>V>U£ - 256U*uf) 

-e 7/2 Q(A^ s + S£b) + Ki 5 - \N x .sBi - ^N^N? - 2U a N£ 5 - - 2N 2 . 5 U* 

+2VjKH + 2K{ 5 V^ - U n K! 2 . - Ki 5 U a - «#V 1J5 - AN^K - ttf + AV 3 Bf b 
-\BUN? + B{) + + MUrNusU* + 8U a U?B l i 5 ) + 0(e 4 ) , (2.18d) 

V 0J = eiUJv - AV™) + e 2 - B^F - K [ 2 ij] + ±6$ - AU a UJ ij 

-WlVt + - e 5 / 2 (±JVi. B< 5« + Xl'l 1 - 

+e 3 Q(iV 2 - B 2 )<^ - 4 Ml + ^ - (JVi^ + U^Ni)** - K 3 2 Ui + AU a K^ ] 

+ 81^%^' - V*{N{ + Si*) + 4iVxyJ^ + AB ik Vj k ' j] 
+AV k ' [i B§ k - AV k B k[i ' j] + lW%UJ ij + 32U a U£Vj - 6AU 2 V^ 

+e 7/2 Q(iV 2 . 5 - B 2 . 5 )5 i3 - Klf + ±J3« - (NuP. + U a N x ^ 

-Kl h Ui + W a Klt + m. 5 V^ + AB ik 5 V} k ^ + AV^B^ + 0(e 4 ) , (2.18e) 

life = e°T) k (U a ) + e 2 Q %(M - BO - 4^ + i(2#' fc + £^ - fl? < 1 )) 

Q °r} fc (iV 2 - £ 2 ) + \{Bf k + B*-* - - N, °T) k {U a ) + - B^ 

-u a ^(ivo + i6f/ 2 °r* fc (c/ CT ) + 4Ki^' fe + 4 °r;. fc (vX) - le^VJ*^ - iw k v^ 

+e 7/2 Q °r* fc (iV 2 . 5 - B 2 . 5 ) + l -{Bf k + Bf£- Bif) 

-JV 1J5 °r} fc (C/ CT ) + Bf. 5 C/^^ - t^) + 0(e 4 ) , (2.18f ) 
where we define 

°T) k (f) = f' k S tJ + f' j 5 ik - f' l S jk . (2.19) 



E. Conversion to the baryon density p* 



We must now convert all potentials from integrals over a, a 1 and o % i to integrals over 
the conserved baryon density p*, defined by Eq. ( 2.10 ). From Eqs. (2/7) and ( p. 12 ), we find 
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a 1 = p*u°v i /y^g , 



a 



p*u°v i v j jy/^g , 



[2.20) 



where u° = (—goo — 2g 0i v i — g^v^)" 1 ! 2 . Substituting the expansions for the metric, Eq. 
( p.4|) , and for the metric potentials Eq. (|2.5|) , we obtain, to the order required for 3.5PN 
equations of motion, 



a = p 



7 4 , 1 .2 



l+e[-v 2 -U a )+e' [-v*+ -v 2 U a - 4v 3 VJ - -N x + -B x + -U* 



4 4 

+2e>' 2 N x ., + e 3 ( + ^v*U a - 10v 2 Vy + 7 -v 2 U 2 + lv 2 N x + 9 -v 2 B x + 



16 



4 



a 1 = pV 



+Au a vy - k& - \n 2 + ~s 2 - w;v; + -u a N x - hj a B x - ^u 3 a 

+e 7/2 (l v 2 N x . 5 + ifl&vV - K[y - \n 2 ., + ~B 2 . 5 - U a N x .^j + 0(e 4 ) 
1 + e ( \v 2 - U„) + e 2 f + lv 2 U a - + -B x - ~iVx + ~£/ 2 



a l ° = p*v % v° 



+2e 5 / 2 iV 1 . 5 + 0(e 3 
1 



<T* = p*V 2 



l + e^-v 2 -U a )+0(e< 
1 



l + e[-v A -U„ 



4 



4 



-2e 5/2 A^ 1 . 5 + 0(e 3 



(2.21a) 

(2.21b) 
(2.21c) 

(2.21d) 



Substituting these formulae into the definitions for U a and the other potentials defined in 
Paper I, Eqs. (4.10) - (4.16), and iterating successively, we convert all such potentials into 
new potentials defined using p*, plus PN corrections. For example, we find that 



U a = U + el-$ x -<!> 2 
7 



+e 2 [^(v*) + -E(v 2 U) - 4£(tW) + -E($x) - S($ 2 ) 

q 1 \ /a (3) \ 

+ 6 -T,{U 2 ) - -S(X)j - e 5 / 2 l-U X" (t)\ + 0(e 3 ) , 



(2.22) 



where henceforth, U, V\ $i, $2, S, and so on, are defined in terms of p* (see Appendix [S]). 



F. Final continuum equations of motion 



Combining Eqs. ( |2.18| ) and ( |2.15| ), substituting the explicit forms of the potentials N x , 
K l 2 , B$, etc. from Paper I, Eqs. (5.2), (5.4), (5.8), (5.10), (6.2), and (6.4), and inserting 
the iterated forms of all potentials, we obtain the equation of motion through 3.5PN order. 

dl! jdt = U' -\- flp^v a 2PN a 2.5PN a 3PN a 3.5PN > (2.23) 
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where 

a PN = v 2 U' { - AvWU* - 3v l U - AUU'* + 8v j V [i ' j] + AV l + \x> 1 + - , (2.24a) 
ai PN = 4vWV^'' fc + v 2 v i U + vV(4$ 2 J ' - 2$f - 2X' j ) - \v 2 {2$i + $f - X'*) 

1 ^ (3) 

+v j v k (2${ k >' - 4$p + 2Pf '* - 4Rp) + ^(3$ 2 - -$! - - X +4V k U' k ) 
+^'(8Vf' j] - 16<Z>f j] + 4X [ij1 + 32G [ j' j] - 16C/y [iJl - 4£' [ >V) + 8V i U> j 

-4$i j C/' j + 8^1/^ + 4F*P + 2£(vV) + 4f/$ 2 * + 4$ 2 [/'* + SP 2 ^ - E>*($ 2 ) 

3 .. .. I .. 3 ... I ... ( 3 ) I W 

+-£-*([/ 2 ) - 2C/X' 1 - 2Xf/' i - 8UV 1 - -£' l pT) + -Xi % - -Xg + 2X 1 +— Y' 1 
2 y 1 2 y ' 4 1 2 2 24 

+4l£ - 8$ 2 - 6G^ - 4G 2 * + 8G'£ + 8G'l - 4Q l + 16G l 7 - 4P?U' j - AH' 1 , (2.24b) 

o (5) i (5) (4) (3) , (3) (3) 

= -x j (l ij — S ij l kk ) + 2v j l ij +2U> j I ij +-U>* X kk -X> ijk l jk 
5 V 3 ; 3 



4<S>[ j - Alf) + -£'> 4 ) + -£'> 2 P) - 4£'V'V°) - -£'*($i) - 6P$f - 2<f> 1 U> 
8 2 2 



l 2.bPN 



2 (5) 2 (4) 

jw + gro j« ? (2.24c) 

15 3 

1 {7) 
4.5PJV = ^(13rV<^' - 4rV<F fc - iVi*) V k 

1 (6) 

+ — (10rV<^ + 4(u • x)x fc <^ - rV5 jfe - 4xVu* + 3x J x k v l ) V k 
30 

1 (5) 

+— (10(u • x)v k 5 ij - x k v 2 5 ij + 2(v ■ x)v l 8 jk + 2a; W* - 5iVu* + 4xVu fc ) V k 

1 (5) 1 (5) 

+— (5r 2 P J - 35^C7 + 59X-n P- 5 + — (5r 2 P'* + 30x l U - 34X'* - 6x J X«) l kk 
15 V ; 30 V ; 

1 (5) 

+ — (15x j x k U' i - Qx j X' ik - 30x l X' jk - 15r 2 X^ k - AY' ijk + 5x l Y' ijkl ) V k 
90 

1 (4) i (4) 

+-v i (v 2 S jk - v j v k ) X jk --(2x j U + 22V j - UX' j - 12v k X> jk ) l ij 

+—(2Ax j v k U' i + 12v i x j U' k + \2x j V Ki + 54v' l X' jk - 72v j X' ik 
18 

(4) i (4) 

+ 12x 3 X' lk + mX 3 ' kl - 72X l ' 3k - 5Y> ljk ) V k — (6v l U - 8V 1 - 2X l ) Z kk 

3 

(3) 

+ (2v 2 U' 1 - 8UU' 1 - 8v k V Kl + 3$f - 2$$ + X' 1 ) T 3 

1 (3) 

+ -(2v 2 U' 1 - 8wVP J - 24PP' 1 - 6v l U + 8V l + 16v J V [lj] + 3$^ - 6$ 2 * + X' 1 ) I kk 
3 

+-(48u J V fc>i - 12v j v k U' 1 - 6v 2 X' ijk + 24v i v l X' jkl - 2AX^ k - §Xt ]k + 6Xi jk 

(3) 

+18u i X J '* - A8v l X [i ' l]jk + 24UX> ijk + 24U' i X> jk - I8${ k ' 1 + Q^(X' jk ) - Y ijk ) V k 
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1 (7) 

+ (10x i x j 6 M - 25x k x l S ij - 9r 2 6 ij 5 kl ) V kl 

630 v ; 

1 (6) 

+— (4x [i v j] 5 kl - 2(v ■ x)5 ij 5 kl - Wx k v l 5 ij ) V kl 

1 (5) i (5) 

(v 2 5 ij 5 kl + 6v i v j S kl + 5v k v l 8 ij ) l jkl +-(4JJ5 3k - 2x j U' k + X' jk ) T jk 

45 9 

1 (5) i (5) 

+— UX' ij - 10x J U 4 ) V kk +—(6x j X' ikl - Y' ijkl ) V kl 
45 V ; 54 v ; 

2 (4) 
+-{U6 ij - t/f/ J - 2v j U' i - - X' ij ) l jkk 

1 (6) 16 (5) 
+— (3r 2 e qik + 2x i x j e qjk + Ax j x k e qij ) J qk x J v k e qjk J qi 

45 45 

2 (5) 
+ — (2(v ■ x)e qik - 2x i v 3 e qjk + ^x j v { e qjk + 12x j v k e qij + Ax k v j e qij ) J qk 

45 



(4) 



--(13Ue qik + x j U' i e qjk - 2x j U' k e qij - X' lj e qjk - 2X' jk e qij - 2x l X' ijk e qlj ) J qk 
9 

(4) 4 (4) 4 (3) 
+ -(Av j v k e qij - v 2 e qik ) J qk -- x j U' k e qjk J qi ~-Ue qtk J qk 

4 (3) i (7) 

+-(v i U' j + 2v j U' i + V j ' { + X' ij )e qjk J qk x* l jjkk 

9 V ; 840 

1 (7) 1 (6) i (5) i (6) i (6) 
j x 3 j-tjkk j i -jjjkk j tj,i q-jjkk j ' qtj jqkk j qik jqjk 

35 40 24 30 15 

1 (5) (5) (5) 1 (5) 1 (5) 

_l vHe q ^ k J qtk —t qtk J q i k —e q% i Ji kk ^ x l J\z[ kk " x^ J\A kk ^ 

15 ^ 30 15 

1 (4) o (4) 1 (3) 2 (3) 
— v * M kkjj +-v j M ijkk +-U A M jjkk +-U' j M ijkk 

6 3 6 3 



i (3) 23 (7) 2 (6) 1 (5) 

—X^ k JvV m — — — — T % ^ kk +— e w Jii kk j\A kk ji l 
3 4200 75 30 



(2.24d) 



Because of their length, we shall defer presentation of the 3PN contributions to later pub- 
lications when they will actually be needed for calculations. 



III. TWO-BODY EQUATIONS OF MOTION TO 2PN ORDER 

A. General treatment of "spherical pointlike" masses 

We must now integrate all potentials that appear in the equation of motion, as well as the 
equation of motion ( |2.24| ) itself over the bodies in the binary system. We treat each body as 



a non-rotating, spherically symmetric fluid ball (as seen in its momentary rest frame), whose 
characteristic size s is much smaller than the orbital separation. We shall discard all terms 
in the resulting equations that are proportional to positive powers of s: these correspond to 
multipolar interactions and their relativistic corrections. The leading Newtonian quadrupole 
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effect is formally of order (s/r) 2 relative to the monopole gravitational potential m/r, but for 
compact objects such as neutron stars or black holes, s ~ m, so effectively this is comparable 
to a 2PN term. Furthermore, if the quadrupole moment is the result of tidal interaction with 
the companion, the size of the induced moment is of order (s/r) 3 , so the net effect is 0(s/r) 5 , 
or roughly 5PN order. Such leading multipolar terms can be calculated straightforwardly, 
but here we ignore them. 

We also discard all terms that are proportional to negative powers of s: these correspond 
to self-energy corrections of PN and higher order. We shall assume that all such corrections 
can be merged uniformly into a suitably renormalized mass for each body, in line with the 
Strong Equivalence Principle. This should be checked by direct calculation, but here we 
ignore such terms. 

We retain only terms that are proportional to s°. For the most part, these are the 
expected terms that depend on the two masses, terms that one would have obtained from 
a "delta-function" approach that discarded all divergent self-energy terms. However, at 
higher PN orders, another class of s° terms is possible, at least in principle. These are terms 
that arise from non-linear combinations of potentials. One could imagine one potential 
being expanded in a multipolar expansion about the center of mass of one of the bodies in 
positive powers of s, multiplied by another potential which is a "self-energy" potential of 
that body, dependent upon negative powers of s. One could then end up with a term that 
has a piece that is independent of the scale size s of the body, but that still depends on its 
internal density distribution. We will show that such terms cannot appear at 1PN order 
by a simple symmetry argument. At 2PN order, terms of this kind could appear in certain 
non-linear potentials, but in fact vanish identically by a subtler symmetry. At 3PN order, 
such s° terms definitely appear, but whether they survive in the final equations of motion 
is an open question at present. We will discuss these matters explicitly at each PN order. 
This approach is, in some sense, a "quick and dirty" version of the Hadamard partie finie 
technique, but with the virtue that finite-size or structure-dependent terms can in principle 
be systematically kept and examined. 

Our assumption that the bodies are non-rotating will imply simply that every element 
of fluid in the body has the same coordinate velocity, so that v l can be pulled outside any 
integral. This assumption can be easily modified in order to deal, for example, with rotating 
bodies. 

Finally we assume that each body is suitably spherical. By this we mean that, in a 
local inertial frame comoving with the body and centered at its baryonic center of mass, 
the baryon density distribution is static and spherically symmetric in the coordinates of 
that frame. In Appendix |B|, we show that the transformation between our global harmonic 
coordinates x l and the spatial coordinates x l of this frame can be written in the form 

x l = x\ + x j {5) + e(A) + B) k x k ) + 0(e 2 )} , (3.1) 

where the subscript A refers to the Ath body and x\ denotes its baryonic center of mass. 
The coefficients A 1 a include the effect of Lorentz boosts, and the coefficients Bj k depend 
on the acceleration of the frame in the field of the companion star. Then, in terms of 
the local coordinates x l , we assume that the baryon density is spherically symmetric and 
static, so that p*(t,x) = ps{r). As a consequence, a finite-size moving body will no longer 
appear spherical, in part because of the Lorentz-FitzGerald contraction. This will result 
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in relativistically induced multipole moments for the body, albeit of order e relative to the 
monopole moment. Ordinarily, these would result in terms of positive powers of s in the 
equations of motion, which we ignore (in other words, as the body's size shrinks to zero, 
the flattening become irrelevant); however, as before, in terms with products of potentials, 
we must worry about the effect of self-potentials with negative powers of s offsetting the 
positive powers from the flattening. We show in the Appendix, however that no such terms 
arise in the equations of motion at 2PN order, but that they will contribute in principle at 
3PN order. 



B. Newtonian and PN terms 

We shall evaluate the acceleration consistently for body #1; the corresponding equation 
for body j^l can be obtained by interchange. At the end, we shall find the centre-of-mass 
and relative equations of motion. 

The Newtonian acceleration is straightforward: 

(a\) N = -(1/mx) j^d'x J P A ^~^} d 3 x' 

= -(1/mi) / / p*p*% ^d 3 xd 3 x' - (1/mO / p*d 3 x / p*'f . (3.2) 

JlJl |x — X | d Jl J2 |X — X'| d 

The first term vanishes by symmetry, irrespective of any relativistic flattening or any other 
effect (Newton's third law). Substituting Eq. ( jSTlf ) for each body and expanding the second 
term in powers of x, using the general formula 

r^ = E^V«(-) , |y|<|x|=r, (3.3) 



<?=0 



we find that all contributions apart from the leading term are of positive powers in s, 
including the effects of relativistic flattening, and thus are dropped, with the result 

(a\) N = -m 2 nyr 2 , (3.4) 

where we define x = xi — x 2 , r = |x|, n = x/r. 

The 1PN terms are similarly straightforward. A term such as v 2 U ,% is integrated over 
body 1 by setting v = V\ and writing U — U\ + U 2 . With v 2 pulled outside the integral, 
the integration is equivalent to that of the Newtonian term (|3.2j ), with the result v 2 !!' 1 — > 
—m2v\n % jr 2 . Other 1PN terms involving quadratic powers of velocity (v l U, v^V^\ $i and 
the velocity-dependent parts of V 1 and X )% ) are treated similarly. Relativistic flattening 
plays no role through 3PN order. 

In the non-linear term UU' 1 , the term involving U-JJ^ is of order s l /s 4 , where s l represents 
a vector, like (x — x') 1 that resides entirely within the body. In Appendix [TJ we argue that 
relativistic flattening introduces corrections of order e(A+Bs k ), e 2 (C+Ds k +Es h s l ), with the 
generic term scaling as e n (l + s) n . The leading term (s*/s 4 ) vanishes by spherical symmetry. 
Any contributions of overall order s~ 2 , s _1 or s +n are discarded. The only way to get a 
term of order s° out of s*/s 4 is to have a correction term of order s 3 , which is automatically 
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of order e 3 , which results in a 4PN term. In the two cross terms U1U2 and U 2 Ui, U\ 
and U( are of order 1/s and s l /s 3 respectively; expanding U 2 about the center of mass of 
body 1 using Eq. ( |3.3| ) yields only products of vectors s Q , including the contributions from 
relativistic flattening. Thus the only terms in the product that vary overall as s° will have 
odd numbers of vectors s l , whose integral over body #1 vanishes by spherical symmetry. 
Only the term from U 2 U 2 contributes, and relativistic flattening produces only corrections 
of positive powers of s. The result is UU' 1 — > —m\n l /r z . 

In the terms V % and X' 1 , the acceleration dv l /dt appears. Working to 1PN order, we 
must insert the Newtonian equation of motion; but working to 2PN order (or higher), we 
must insert the 1PN (or higher) equations of motion; the 2PN terms so generated will be 
discussed in the next subsection. For V 1 , the result using the Newtonian equation of motion 
is 



V 



p*"{x' 



X 



d 3 x'd 3 x" + 



X" 



p*W • (x — x') j3 



,/|3 



dV. 



(3.5) 



The double integral is integrated over body #1 similarly to the term UU' 1 , and the velocity- 
dependent term is integrated similarly to the term v 2 U' 1 . The general result of these consid- 
erations is that, at 1PN order, only terms are kept in which, in the quantity x — x', the two 
vectors are evaluated at the baryonic center of mass of the two different bodies, respectively, 
and never within the same body. 

The resulting N and 1PN equation of motion is 



a 



1 (PN) 



m 2 i , m 2 i 



a 



2(PN) 



. m 2 , 
#2. 



I — + 5 v{ + 4 Vl ■ v 2 - 2v 2 2 + -O2 ■ n) 

r r 2 



t> 2 )*(4f2 • n — 3t>i • n) 



(3.6) 



Under the interchange 1 ^ 2, n l 



-n 



C. 2PN terms 

Since we are only working to 2PN, 2.5PN and 3.5PN orders here, we may evaluate the 
2PN terms without regard to the effects of relativistic flattening, since the corrections would 
be of 3PN or 4PN order or higher. We only need to expand any potentials about the centers 
of mass of the bodies to identify all possible s° terms. The terms explicitly cubic in v 1 [the 
first two terms in Eq. ( |2.24| b)1 are simplest to evaluate, since the potentials involved already 
appeared at 1PN order. Integrating over body #1, discarding terms that are singular in s, we 
obtain a* 1 (l&2) = — (m 2 /r 2 )[4(t>i -t^Xfi -n) — vf(v 2 ■ n)]v{. The next ten terms are explicitly 
quadratic in velocities; the potentials are both linear and quadratic in the masses, and some 
involve time derivatives that require substituting the Newtonian equation of motion. Pulling 
the velocities outside the integrals leaves potentials similar to those that have already been 
integrated at 1PN order. The potential P%* also appears; unlike the potentials encountered 
so far (U, V\ X, $ 2 , . . .), which depend on the pairwise separation between points, i.e. 
on the distance tab = |xa — x#|, P$ and others like it depend on the distances between 
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points taken as a triplet, namely on the combination A (ABC) = tab + ^ac + r BC- We 
denote P% and similar potentials like Gi, G2 etc. as "triangle" potentials. Their evaluation 
is discussed in Appendix |Cj. From these terms only the obvious "point" mass terms arise, 
while all others are proportional to either negative or positive powers of s. We obtain, for 
example, v 2 Q 2 l ~~ > —mim2V 2 ri l / r 3 , v^v k ^' k — > —m<iv\{y\ ■ v<i){y\ ■ n)/r 2 and v^v k P^' k — > 
mim2[n l (8(vi ■ n) 2 — 3vf) — bv\(vi ■ n)]/4r 3 — m|[n*(4(i;i ■ n) 2 — vf) — 3v\(vi ■ n)]/Ar 3 . 

The next 13 terms in Eq. ( |2.24j b), which are explicitly linear in velocities, are quite 
similar, except that they involve either an additional time derivative of potentials, or vector- 
like potentials (V 2 \ $ 2 5 Triangle potentials appear in several places (G 7 , P^)- As before, 
expansion about the centers of mass yields the normal point mass terms and terms of only 
positive and negative powers of s. 

Of the remaining 33 terms, many involve several implicit powers of velocity coupled to 
potential-type expressions; integration of these terms over body #1 is handled as before. 
However, several terms do not involve velocity at all, and thus are cubically non-linear in 
masses. Examples include the terms U 2 U' 1 , S' i ($ 2 ), and the portions of the terms XU , ' t , 
Xi, y- 0000 ^ where time derivatives have generated accelerations, thence the Newtonian 
potential. In these cases, the high-degree of nonlinearity presents at least the possibility of 
new contributions at order s°. The term U 2 U ,% illustrates the issue. Writing U — U\ + U2, 
we have (Uf + 2L r i£/ 2 + U 2 ){U( + Utf), to be integrated over body #1. Consider the term 
Uf U2 ■ Expanding U 2 about the center of mass of body #1 using Eq. (|3.3j ) and integrating 
over body #1, we obtain 

mi A k |x-x'| Ji |x-x"| Y 9! Vr/ v ' 

where the "barred" coordinates are all defined relative to the center of mass of body #1, 
and r = |xj. — X2I. Only the term with q = 2 can produce a contribution of overall order s°. 
Using the spherical symmetry of body #1, the result is 

i_L Qf p *f 2 U 2 d 3 x^j V k V ijk ^) oc V 1 V 2 ( J) = ■ (3-8) 

Note that the integral in Eq. ( |3.8j ) scales as s° for a fixed mi, yet depends on the internal 
structure of the body. Nevertheless, the term vanishes via a combination of the spherial 
symmetry of the body, and the fact that V 2 (l/r) = 0. Other possible s° terms also vanish 
by symmetry, with the final result that U 2 U' 1 — > —m\ri l /r 4L . Similarly, for example E'*($ 2 ) - > 
—mim\n l / V 4 . 

Additional cubically nonlinear terms are P^U'i, and the acceleration-generated terms 
in G 2 l , G'l and G-j 1 . Since these involve the triangle potentials, they will be discussed 
in Appendix |Cj they generate no s° terms, again because of symmetry combined with 
V 2 (l/r) = 0. The term H>< = V' i P(C/' J ' fc P| fe ) involves a still more complicated "quadrangle" 
potential, which is a function of four points. In Appendix |C| we show that it likewise 
generates no s° terms, with the result 

(1/mi) jf p*H*d 3 x = ^ [2m im2 + ^ ) . (3.9) 
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Working to 2PN order, we must also include terms generated by substituting the 1PN 
equations of motion into accelerations generated by time derivatives acting on velocities in 
1PN potentials, specifically the terms 4V l + X ,l /2. This leads to the integral 



— / p*d 3 x 
rrii Ji 



^(^ + ^>yi,x') 



(3.10) 



where h % = (x— x') 1 / |x— x'|, and we substitute Eq. Q2.24|a ) for a 3 PN (t, x'). Evaluation of these 
terms follows the same method already outlined for normal, non-triangle, 2PN potentials. 
The resulting 2PN contributions to the one-body equation of motion are 



a l (2PN) 



a 2(2PN) 



m 2 



-n 



— ~ {^ v 2 ~ 8^1 ' v 2 + 2(t>i ■ nf — 4(t>i ■ n)(v2 ■ n) — 6(t>2 • n)' 



15 39 

-[ ~ v 2 ~ t; v i ' v 2 — t v i + ~2^ Vl ' ~ 39(^1 ' n)(v 2 ■ n) 



mi ( 5 2 5 

57 m 2 69niini2 

4 r 2 2 r 2 
3 



17, , 2 
— {v 2 -n) 



m 



2v\ + Av 2 2 (vi ■ v 2 ) - 2(vi ■ v 2 y 



+^v 2 (v 2 ■ n) 2 - 6(ui • v 2 )(v 2 ■ n) 2 + ~v%(v 2 ■ n) 2 - ^-(v 2 ■ nf 

2 2 o 



, m 2/ i 
+ 72-K 



nil 
4r 



^55(^2 ■ n) - 63(f 1 • n)^j - (^(vi ■ n) + (v 2 ■ n) 



+v 2 {v 2 ■ n) + Av 2 {vi ■ n) — 5v 2 (v 2 ■ n) — 6(v x ■ n)(v 2 ■ nf 

9 

—A{v x • n){yi ■ v 2 ) + A{v 2 ■ n)(v x ■ v 2 ) + ~(v 2 ■ nf 
1^2. 



These agree completely with other results 



(3.11) 



IV. RADIATION REACTION TO 2.5PN AND 3.5PN ORDER 



Since the multipole moments are strictly functions of time, integrating the 2.5PN and 
3.5PN terms in Eq. ( |2.24| ) over body #1 is straightforward. The 2.5PN terms are either 
trivial, involving x l or v l , or are similar to integrating the Newtonian term. The result is: 



a l (2.5PN) 



a 2 (2.5PN) 



o . (5) 1 

-x{(V j — 5 13 T 
5 U 3 

(5) 



(5) (4) , (3) (3) 

^ + 2v{ l ij ---^-n 1 X kk -3^^n% fc V k 



3 r 1 



2 {0> 2 ^ 
15 + 3 C J ' 



(4.1) 



Likewise, the 3.5PN terms are either trivial, or involve integrating Newtonian or lPN-like 
potentials. To keep the expression for the 3.5PN terms simple, we assume that, to lowest 
order, miXi + m 2 x 2 = (see Sec. |V A| for discussion), so that we can write, in the 3.5PN 
term only, x\ = (m 2 /m)x t and x % 2 = —{nii/ni)x l . Defining a 2 = ni2/m, we obtain 
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1 (7) 

al(3.5PA0 = ^a 3 2 r 3 (13n k ^ - 4nW k - n^n k ) V k 

1 (6) 

+— a 3 2 r 2 (Wv k 5 ij + 4rn k 5 i] - v l 5 jk - 4nW + 3n j n k v l ) V k 
30 

1 (5) 

+— a 3 r(10fv k 5 ij - n k v 2 5 lj + 2rv i 5 ik + 2n l v 2 5 jk - 5n i v j v k + An j v l v k ) V k 



(5) 

+— to 2 {(55 - 3Qa 2 )n k 5 ij - (19 - 10a 2 )n i 8 ik + (2 + Qa 2 - 10a 2 2 )n i n j n k \ 2 jk 
15 

1 (4) 



+-a 3 v i (v 2 5 jk -v j v k ) l jk 
1 m 2 



{ 42 (1 _ a 2 )v k 5 lj - 6(3 + a 2 )rn fc <^' + (7 - a 2 y5 jk + (1 - a 2 )rW fc 
6 r L 



(4) 



-(19 + 3a 2 )v i n j n k + 2(15 - 3a 2 - 4a 2 Vn i n fe - 3(5 - 9a 2 + Aa 2 2 )fn i n j n k } V k 
+^TT l 16 ^ + «2)-nVn fe + (2 - c^n^u* - 6wVn fc ) 

+ (3 + « 2 )n fe (6rwV - 4uV) + 3(1 - « 2 )Vn fc (5rV - 4tV)} J jfc 

--^ (8(1 - c^-ra* - 2(2 - a 2 >V + 2(3 + a 2 )rv l + 3(1 - a 2 ) Wl J' 
6 H I r " J 



(3) 



1 



(7) 



+ a 2 2 r 2 mn i n j 5 kl - 25n k n l 8 ij - 95 ij 5 kl ) X jkl 

630 2 V ; 

1 (6) 

+— a 2 2 r(4n [i v j] 5 kl - 2r5 ij 5 kl - 10n k v l 5 ij ) V kl 

1 (5) 

a 2 2 (v 2 5 ij 5 kl + 6v i v j 5 kl + 5v k v l 5 ij ) V kl 

45 

1 (5) 

+ — r ^l {(43^ + 5n k n l )5 ij + (7 + Wa 2 )n i n j 5 kl - 15(1 - 2a 2 )n i n j n k n 1 } 
90 r 1 J 



2 TO 

+ £ (t)V + 2a 2 v j n i - 3(1 - a 2 )rn i n i \ V 

1 (6) 16 (5) 
+— a\r 2 (3e qik + 2n i n j e qjk + An j n k e qij ) J qk —f-a 2 rn j v k e qjk J qi 

45 45 

2 (5) 
+— a 2 2 r{2fe qik - 2n i v j e qjk + 5n j v i e qjk + 12n j v k e qij + An k v j e qij ) J qk 

45 

2 (4) 2 to , , (4) 

+-a 2 2 Uv j v k e qij - v 2 e qik ) J qk --— {(1 + a 2 )n l n j e qjk + 10e qik + 2n j n k e qi H J qk 
9 9 r 1 > 



[v l n> + 2a 2 v^n l - 3(1 - a 2 )rn t n : >} e qjk J qk 



4TO2 • 

'972" 

1 (7) 1 {?) 1 (6) 1 TOo (5) 

— a 2 rn { l" kk +—a 2 rn j V lkk +—a 2 v i X" kk -^r-^n* I jjkk 
840 35 40 24 r 2 
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1 (6) i (6) •, (5) (5) 

a 2 rn 3 e qij J qkk a 2 rn 3 e qik J qjk a 2 v j (2e q[llk J^ j]k +e qi] J qkk ) 

30 15 15 

1 (5) i (5) i (4) 2 (4) 

-— aW M kkjj -—a 2 rn J M kklJ --a 2 v l M kk " +-a 2 v 3 M ijkk 
30 15 6 3 

1 m (3) (3) 90 (7) 2 (6) -, (5) 

+ — l n HM jjkk -6n j n k M jkU ) l ijjkk +—e qij J qjkk +— M kkjj \ 

6 r 2 v ; 4200 75 30 

a 2(3.5PN) = 1^2. (4.2) 



V. RELATIVE EQUATIONS OF MOTION 
A. System center of mass and the transformation to relative coordinates 



It is useful to note that the 1PN equations ( |3.6|) , including the 2.5PN terms (}4.1|) , admit 
a first integral that corresponds to uniform motion of a "center of mass" quantity, namely 

mi vi(l + -vl) + m 2 v i 2 (l + -v 2 2 ) - -^^{yi + v 3 2 )(5 ij + n l n 3 ) - mV = & . (5.1) 



where 



= _ e 9*fc 

15 3 



(5.2) 



C % is a constant, and where we have assumed that, to Newtonian order, m\x\ + m 2 x 2 = 0. 
The 2PN corrections to this first integral will not be needed here. Equation (|5.1| ) can also 
be obtained by calculating the system dipole moment T 1 to the corresponding order (see 
Appendix p 2| ). Choosing the coordinates so that C l = 0, we obtain the transformation 
from individual to relative coordinates and velocities, to 1PN and 2.5PN order, 



i m 2 i 1 5m. 2 m ■ ■ 

i = — x + -1] — (v )x + V , 

m 2 m r 

i m i i 1 5m, ~ m. , 

x\ = -x l + -7] — (v 2 )x l + V 1 

m 2 m r 

j m 2 ■ 1 5m 

V\ = V 1 + -71 

m 2 m 

4 m l 4 1 $ m \, 2 m \ 4 

v 2 = — v i + -ri — (v 2 — y 

m 2 m [ r 



m 



-rx 



m 



-rx 



(5.3) 



where 5m = m% — m 2 . These transformations do not affect the Newtonian term, of course. 
However, the 1PN and 2.5PN corrections in Eqs. (|5.3j ) will generate 2PN and 3.5PN terms, 
respectively, when we transform the 1PN terms in the equation of motion to relative coor- 
dinates. The multipole moments that appear in the 2.5PN terms in the equation of motion 
Qi.lD must also be converted to relative coordinates, keeping any PN corrections generated 
by Eqs. (|5.3|); this is treated in Appendix p 1|. In addition, in the 2.5PN terms, multiple 
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time derivatives of the multipole moments will generate accelerations, for which the 1PN rel- 
ative equations of motion must be substituted; in the explicitly 3.5PN terms, the Newtonian 
equation of motion suffices. 

Calculating a\ - 4 using Eqs. (U), flPD and (gj), substituting Eqs. and 

the time-derivatives of the multipole moments (|132| ) , we obtain the final relative equation of 
motion through 2.5PN order, plus 3.5PN terms as given in Eqs. ( |1.2|) and (|1.3| ). 



B. 3.5PN radiation reaction and energy-angular-momentum balance 

A useful check of our radiation-reaction terms at 2.5PN and 3.5PN order is to verify 
that the resulting energy and angular momentum loss in the orbital motion is identical 
to previously derived energy and angular momentum flux expressions, accurate to 1PN 
order beyond the quadrupole approximation. In fact, Iyer and Will |10|JTT[| approached this 



from the opposite direction, beginning with the 1PN accurate flux expressions [EB-pSj, and 
deriving the most general form of a two-body relative equation of motion at 2.5PN and 3.5PN 
order required by energy and angular momentum balance. Writing the radiation-reaction 
terms in the equation of motion ( |1.2|) in the general form 

M.bPN = a-iv 2 + a 2 m/r + a 3 r 2 , 
B 2 .5pn = biv 2 + b 2 m/r + b 3 r 2 , 

A3.5PN = Cif 4 + c 2 v 2 m/r + c 3 v 2 r 2 + c^r 2 mjr + c 5 f 4 + c 6 (m/r) 2 , 

B3.5PN = d\v A + d2V 2 m/r + d 3 v 2 r 2 + d^mjr + d^r 4 + d§{mjr) 2 , (5.4) 

they showed that energy and angular momentum balance would hold if and only if the 
coefficients a^, bi, Ci and di satisfied the following equations: 

a\ = 3 + 3/5 , a 2 = 23/3 + 2a — 3(3 , a 3 = —5/3 , 

b\ = 2 + a , b 2 = 2 — a , 63 = — 3 — 3a , (5.5) 

for the 2.5PN coefficients, and 

c a = -L(117 + 1327/) - ^/3(1 - 3r7) + 36 2 - 36 6 , 

c 2 = -^(297 - 310r/) - 3a(l - 4r/) - |/3(7 + I3rfj - 28 x - 35 2 + 35 5 + 35 6 , 

c 3 = ^(19 - 7277) + ^(1 - 377) - 55 2 + 55 4 + 55 6 , 

c 4 = - — (687 - 36877) - 6ar] + -/3(54 + 17r/) - 25 2 - 55 4 - 65 5 , 
28 2 

c 5 = -75 4 , 

c 6 = -^(1533 + 498t/) - a(14 + 9r/) + 3/3(7 + 4t/) - 25 3 - 35 5 , 

d 1 = ~(2 + a)(l-3r ] )-5 1 , 

d 2 = -4(!39 + 768r/) - i a (5 + 17r/) + 61 - S 3 , 
o4 Z 
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d 3 
d i 
d 5 



^(369 - 62477) + j|(3a + 2/3)(l - 377) + 3^ - 3<5 6 , 

^(295 - 33577) + i«(38 - II77) - 3/3(1 - 377) + 2Si + 45 3 + 35 6 . 
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(19 - 7277) - 5/3(1 - 377) + 55 6 



• — (634 - 6677) + a{7 + 377) + 63 



(5.6) 



for the 3.5PN coefficients. The two degrees of freedom (a, (3) at 2.5PN order and the 
six (5i) at 3.5PN order correspond to gauge or coordinate freedom, and have no physical 
consequences. For example, at 2.5PN order, the values a = —1, (3 = correspond to the 
gauge used by Damour and Deruelle ||, while the values a = 4, (3 = 5 correspond to the 
so-called Burke-Thorne gauge (see for example §36.11 of |5J), also used by Blanchet . 



It is then a non-trivial check of our result to verify that the 18 coefficients in our 2.5PN 

]1.3c[) and (|1.3d| ) yield a unique, self-consistent solution for the 8 



and 3.5PN terms, Eqs. 
gauge coefficients. The result is 

a= -1, (3 = 0, 

— +677, 

10 , 5 5 = 



77 

do = -— - 



5 242 

1 n 

42 21 ' 



3 


5s = 


79 


92 


2^' 


14 ~ 


' y 


5 6 = 


439 


18 




28 


+ T 


77. 



(5.7) 



VI. CONCLUDING REMARKS 

We have successfully used DIRE to derive equations of motion for compact binary systems 
through 2.5PN and to 3.5PN order, with results consistent with other methods. Instead of 
using formal delta-function, matching, or regularization techniques to treat the bodies, we 
modeled them as fluid balls, considered to be suitably spherical and non-rotating, and small 
compared to their separation, and carried out explicit integrations over them. We then used 
a technique whereby we could identify those contributions to the equation of motion that 
are independent of the scale size of the bodies (for given masses). This method can be 
extended straightforwardly to the complicated 3PN contributions, to spinning bodies, and 
to bodies with tidal interactions. We can also consider the effects at higher PN order of 
internal self-gravity (contributions of order s~ n ). These are the subjects of ongoing research. 
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APPENDIX A: KEY FORMULAE USED IN THE EQUATIONS OF MOTION 



Here we summarize some of the key formulae from Paper I |2| that will be needed here. 
The potentials that appear in the equations of motion are all Poisson-like potentials and 
their generalizations, namely a superpotential and a superduperpotential: 

fit,*') 

IM |x 



P(f) = ^f l^-^x', V 2 P(f) = -f, 
47T Jm |x — x'l 

S(f) = ^f /(t,x')|x-xVV, V 2 S(f) = 2P(f) 
Air Jm 



SD(f) = ±- ( /(i,x')|x-x'| 3 dV, V 2 SD{f) = l2S{f). (Al) 

47T JM 

Note that, in evaluating Poisson potentials and superpotentials of sources that do not have 
compact support, our rule is to evaluate them on the finite, constant time hypersurface Ai, 
and to discard all terms that depend on the radius of the near-zone, TZ. Unlike Paper I, we 
now define all potentials in terms of the conserved baryon density p*: 

Jm |x — x'| 

X(f)= [ p*(t,x')/(t,x')|x-x / |rfV = 5(4 7 r P 7), 
Jm 

Y(f)= [ p*(t,x')/(t,x')|x-x'|W = 5 J D(47r P 7)- (A2) 
Jm 

The specific potentials used in the 2PN, 2.5PN and 3.5PN equations of motion are then 
given by 

V i = E(u*) , $1 J = E(«V) , 
$ 2 = £([/), X = X(1), 

$*ee£(U*), Y = Y(1), 
X, = X(v 2 ) , X 2 = X(C/) , 
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$i 


= £(„ 2 ), 


U' 

v 2 






EE , 




EE P(U^U J ) , 


G l 


s P(f7 2 ) , 


G 3 


EE -P((j> k V k 


G 5 


EE -P(V k U' k 


G\ 


EE P(f/' fc U fc ' i ) 



P, = p« = $ 2 _ -U 2 



2 

G 2 EE P(C/[/) , 

G 4 ee p(y^V^) , 
G 6 ee P(£7«$?) , 

- -P{U^U) , 77 ee P{U' ij Pi j ) . (A3) 

The multipole moments that appear in 2.5PN and 3.5PN terms are defined by 



P^= [ 


r^d 3 x , 


Jm 




J Q EE / 


T°°X Q d 3 X 


Jm 




rjiQ ^iab 


[ T° b X aC 




Jm 


M ijQ = f 


r ij x Q d 3 x 


Jm 





(A4) 
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To the order needed for our purposes, 

T 0O = a _ a H + 4<TUff _l_ VU 2 

r°* = a< + 4*7*17, + -U*V™ + -^itf , 

7T 47T 

r« = <r« + -L (ttfttf - ^ Vf/ ff 2 ) . (A5) 



APPENDIX B: TREATMENT OF "SPHERICAL" BODIES IN PN EXPANSIONS 

We define our bodies to be spherical in a suitably chosen comoving frame. For a given 
body A, we choose a frame that momentarily has the same coordinate velocity v relative to 
the global PN frame as body A. Also, in the limit tua — > 0, the frame is locally Lorentzian 
with its origin at x^, i.e. the frame is a local, freely falling frame in the field of the other 
body. In that frame, with coordinates (t, x l ), the conserved baryon density distribution of 
body A is taken to be spherically symmetric and static, i.e. p*(t,x l ) = ps(|x|). We define 
the baryonic mass, center of mass and velocity of body A according to Eqs. ( |2.13 ). 

Our goal is to calculate PN potentials and to integrate them over one of the bodies using 
the fact that p* is spherical and static in the local "hatted" coordinates. The general form 
of the integrals to be evaluated is / J p*(t, x)p*(7, x')/(x, yi')d 3 xd 3 x'. First we note that the 

quantity p*(t,x)d 3 x = p(t, ~x)u^J (—g)d 3 H fl is a scalar, i.e. is the same at a given event in 
any coordinate system. Thus we only need a transformation of the integration variables x 
and x' in /(x, x') to the hatted coordinates. Notice that x and x' are taken at the same 
global coordinate time t, but will not necessarily be at the same local coordinate time t. 

Consider the event V inside the fluid at x' and the center of mass event Q at xa, both 
at the same time t as the field point x (see Figure 1). The field point could be within body 
A itself, in free space, or in the other body. Define 3 e~, where e„ are the 

basis vectors of our global PN frame; x is purely spatial in this frame. In the local comoving 
frame 

x = x J ej = x^e^ = teg + ffe-- . (Bl) 

Let the transformation between the two basis vectors take the general form 

e, = (Aj + B%x v )e & , (B2) 

where and is symmetric on the lower indices. The coefficients A^ 

correspond to boosts and coordinate re-scalings, while the coefficients correspond to 
making the frame freely falling. Substituting Eq. (|B2|) into (|B1| ) we obtain 

i = \ i y + B%x j x k , 

x i = A{ x k + Bi lX k x l . (B3) 

Notice that V is not simultaneous with Q in hatted coordinates, and that the time difference 
t depends on xP . However, since we assume that ,t) is time-independent, this variation 
of t with integration point is irrelevant. 
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Writing A^ = + A", it is straightforward to show that, in a PN expansion, the 
various coefficients in Eq. (B2) have the leading orders A° ~ A k ~ e, Aq ~ A° ~ e 1 / 2 , 



5q- ~ i?oo ~ B 3 kl ~ e, 5^ ~ ~ 5q ~ e 3 / 2 . The transformation then takes the 

form & = xi + A k x k + B kl x k x l , where the corrections are of leading order e. Inverting the 
transformation iteratively, we obtain 



x 



x\ + V{5) + e(A) + B) k x k ) + 6 2 (q + D] k x k + E) kl x k x l ) + 0(e 3 )} , (B4) 



where the coefficients are functions of the A^ and B kl . To the PN order at which we are 
working, their explicit forms are not needed. Notice that, in terms of the scale s ~ \x\, 
the flattening correction terms in Eq. ([54]) have the general form e n (a + (3s + . . . + 7s") ~ 
e n (l + s) n . 

In the double integral of /(x, x'), where the dependence is generally on the difference 
x — x', there are two cases to consider, one where x and x' are in different bodies, and the 
other where they are in the same body. In the former case, we substitute Eq. ( |B4| ) for 
both x* and x', and expand about the quantity — xg in powers of s using Eq. ( |3.3| ). 
convert the quantities p*d 3 x to the hatted coordinates, and then integrate over the spherical 
density distributions, keeping only terms of 0(s°). The first term in the expansion of Eq. 
( |B4|) produces the normal multipole expansion of the potential, while the remain terms are 
relativistic flattening corrections. In the case where both x and x' are in the same body, we 
have 

(x - x'Y = (x- x'Y{5) + e\A) + B) k [x + x'f\ 

+e 2 [q + D) k [x + x') k + E) kl (x k x l + x k x 1 ' + x k 'x 1 ')} + 0(e 3 )} . (B5) 

In this case, the corrections come only from relativistic flattening, and also have the general 
form e n (l + s) n . 

APPENDIX C: EVALUATION OF NONLINEAR 2PN POTENTIALS FOR 

TWO-BODY SYSTEMS 

1. Triangle potentials 

The potential P% = P(6 r,l 6 rj ) represents a new kind of potential that first appears at 
2PN order. Unlike the Newtonian potential U, whose fundamental ingredient l/|x — x'| 
depends on the field point and the source point, the fundamental ingredient of depends 
on the field point and on two source points, hence the name "triangle" potential. The 
potentials Gi, i — 1 . . . 6 and G\ are also triangle potentials. To see this, we write P^ in the 
form 

1 r ri 3 r' r (i J — T l '\ i r (V — <r w V 

p» = 1. / -—— / P *"d 3 x"f — / p*"'d 3 x mK - — 



4-7T Jm |x — x'| J ' |x' — x"| 3 J |x' — x'"| 3 

E I PA^ A d 3 x A f p* B V B d 3 x B ±- I ^ 

Ti^JA JB 47T Jm x — x' — X' Xb 



A.B 



E /„ ip*AP B d?x A d?x B V l A V B G{xAB), (CI) 

A.B ' 



A JB 
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where 



g(xAB) = - In A(xAB) + 1, 

A(xAB) = |x - x A | + |x - x B \ + \x A - x B | . (C2) 



The triangle function Q(ABC) satisfies the equations 
V 2 A g(ABC) = -l/(r AB r AC ), 
V A V l B g{ABC) = 



1 




( l - 


- 1 ) 


1 1 ■ 


2 


r AB 


\r AC 


r BC J 


r AC r BC . 



(C3) 

together with the obvious results obtained by interchange of indices. Specific gradients of 
Q(xAB) have the form 

Vi5MB) = (Vb + n AB )> , 

V' A V B S(xAB) = A , xAB)2 (VA - n AB Y(y B + n AB )' + —1-—^ - n\ B n s AB ) , 

V b^ A B) = ^f+m + n AB f(y B + ^ - ^^ (^ " -A^b) 

-^IS)^-^' (C4) 

where = (x — = \ya\, and y\ = y\/y A . With these definitions, the other triangle 
potentials may be written 



G 1 


= E 

A,B 






(C5a) 


G 2 


= E 

A,B 


/, 


[ B P A p B d 3 x A d 3 x B (a* B V B + v B v j B V B V j B )g{xAB) , 


(C5b) 


G 3 


= E 

A,B 


/, 


^p\p B d z x A £x B v\v B V B W B g{xAB) , 


(C5c) 


G A 


= E 


/, 


f B P A p B d 3 x A d 3 x B v A viW j A W B g(xAB) , 


(C5d) 


G 5 


= E 

A,B 


/, 


j B P A p B d?x A d*x B {a\V B + v A v A V j A V B )g(xAB) , 


(C5e) 


G G 


= E 

A,B 




f B P A p B d 3 x A d 3 x B v A v A V B V B g(xAB) , 


(C5f) 


G\ 


= E 

4 R 




f B P A p B d 3 x A d 3 x B (v A - h B )V A V B g(xAB) . 


(C5g) 



2. The triangle potential 

For a two-body system, we may write = ^2(11) +^-^2(12) + -^2(22)' wnere the subscripts 
in parentheses denote the contributions from the bodies. The contributions -P^n) an d ^2(22) 
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from a single body with a spherically symmetric mass distribution p* = p(r) can be derived 
in a simple manner. For spherical symmetry, the equation for P 2 takes the form V 2 P 2 = 
—n ll \U') 2 ^ where U' = dll/dr = —m(r)/r 2 , with dm(r)/dr = A7cp(r)r 2 . Writing 

P\ 3 = (l/3)5 lj P 2 + n <l3> PT 2 , (C6) 

where n <y> = n l] — 8^/3, it is easy to show that the trace P 2 and the traceless part PT 2 
satisfy the equations V 2 P 2 = -(U') 2 , and V 2 PT 2 - 6PT 2 /r 2 = -{U') 2 . Demanding that 
the solutions be regular at the origin and vanish at infinity yields 



P 2 = - 1 / p(r')m(r')r'dr' + 4tt p(r')m(r')dr' , 

2H r Jo Jr 

PT 2 = - /„ P( r ) m ( r )r dr + — I — dr , (C7) 



4r 2 5r 3 jo 5 Jr r 

where R is the radius of the body Inside the body, P 2 ~ (m/s) 2 . Outside the body, 
P^^ — (m/2r) 2 (n tJ — 5 lJ ), neglecting internal-structure terms of 0(s _1 ) and 0(s +1 ). 

For the term P 2 (i2)) an< ^ ^ or a P°i n t between the two bodies, we combine the definition 
( |C1|) with the appropriate gradient of Q from Eq. ( |C4|) , and show that only point mass terms 
contribute. Then, for an interbody field point, we have 

1m x m 2 ( {y x -n) {i (y 2 + n) j) 5 ij - n ij \ 
+ A(xl2) ^ A(xl2) r ) ' 1 j 



where n % = x\ 2 /r and A(xl2) = yi + y 2 + r, and where self-energy terms of 0(s x ) and 

2 



0(s +1 ) have been dropped from P| 



3. Evaluation of triangle terms for two-body systems 

We then take either spatial derivatives (e.g. G'l) or time derivatives (e.g. G\) of the 
triangle potentials, and in some cases multiply them by other factors (e.g. v 3 v k P 2 ,,k or 
U' 3 P 2 ), and integrate them over the density of body #1. Because we are dealing with a 
two-body system, then either two or all three of the points x, A and B in Q(xAB) will 
reside in the same body, in other words, in Eqs. ( |U5|) , we encounter the possibilities £?(111), 
£(112), £(121), and £(122). The quantity £(111) and its derivatives are purely internal to 
one body and can be treated fairly simply. Notice that a single derivative of £(111) is of the 
general form s J /s 2 , while a double derivative is of the form s lJ /s 4 , and a general derivative 
is of the form s®/s 2q , where = s ll - tq . Since we must be left with one spatial index i, 
such purely internal terms have odd parity, and must integrate to zero. For the mixed-body 
cases, we must expand the functions £ about the centers of mass of each body, sort the 
terms in powers of the scale s for each body, and retain only final contributions of order s°. 
This will be aided by a general expansion of the function Q(ABC) in powers of tab^ac-, 
where points A and B are assumed to lie inside one body, and point C is inside the other 
body, so that r^B ~ s ^ t ac- Straightforward methods lead to the expansion 
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1 00 ( —r n \ m+1 

G(ABC) = -\nr AC + 1 -l n 2 + -£ 1 [h 



2^o {rn+l)\ 

,{ M ^ A (JL) + ^.rvr (^)} . (09) 

Note that each term in the expansion is of order r AB ~ ~ s m+1 , and depends on gradients of 
1/tac- Since A and C are in different bodies, = 0; this fact will be important in the 

considerations to follow. The first term in the braces produces a contribution of order s m+l 
but of parity (number of n AB vectors ) (— l) m ; no matter how many gradients are taken 
with respect to any of the variables, this relationship will be unchanged. Hence a term of 
order s° will have odd parity, a term of order s _1 or s +l will have even parity, and so on. 
(Because of the additional scalar factor r A B, enough gradients with respect to A or B will 
generate terms of negative powers in s.) The second term in the braces has parity (— l) m+1 
and order s m+l , a relationship again preserved under any gradients. Furthermore, gradients 
of this term will yield terms either completely independent of s or of positive powers in s; 
no negative powers of s or terms proportional only to the unit vector be produced 

by this term. 

Armed with these characteristics of Q, consider as an example, the term G'l in the 
equation of motion. Taking the gradient of Eq. (p5]a) with respect to x, pulling out the 
velocities, integrating over body #1, and considering all the possible cases for A and B, we 
first find no contribution from V^.V^V^(111) ~ s t s : >s k /s 6 , by symmetry. Considering the 
other cases, £7(112), £7(121), and (7(122), we find that the only s° contributions from the 
gradients of the expansion ( |C9| ) are either odd parity (from the first term ), and thus vanish 
on integrating over body #1 or #2, or are independent of s (from the second term), and 
yield the desired "point mass" result: G'l — > mim 2 [n*(4t;i • nv 2 • n — v\ ■ f 2) — 2v\v2 ■ n — v\v\ ■ 
n]/2r 3 — vn^ln % {y2 ■ n) 2 — ri l v\ — v l 2 V2 • n]/2r 3 . 

Consider as a less straightforward example, the term in G 2 * that depends on the ac- 
celeration. Inserting the Newtonian equation of motion a % B = U B , we must evaluate the 
term 

- / P *d*xY: I I P> B d 3 x A d 3 x B (Ul mt + Ul e jV^ B g(xAB)^ (CIO) 



mi 



where we have split U into a contribution from within body B itself ("hit") and 
from the other body ("ext"). The case A = B = 1 has a purely internal term 
from Ui int V l x V B G(lll), which vanishes by symmetry, and a term Lffg^V^Vs £7(111) ~ 
^Eg^V^VV'Y^s's 3 / 54 ) where we expand the external potential about the cen- 
ter of mass of body #1, and where r = r\2. Only the q = 2 term contributes at 
overall order s°, leading to an integral of the schematic form J 1 (x kl s l / s 4 ) V fc ' J r _1 — > 
(6 kl 5 lj + 5 ki 6 jl + 5 kj 6 u )V klj r- 1 oc V 2 VV -1 = 0. Similarly, combining u£ tira ~ s j /s 3 with the 
expansion of V l x V 3 B Q(112) produces a potentially s° term only from an order s 2 term from 
the derivatives of Q; but such terms are necessarily accompanied by several (three or more) 
gradients of r _1 ; because only a single index i remains at the end, two of the gradients are 
always contracted into V 2 , which vanishes when acting on r _1 . All the possible combinations 
of A and B in Eq. (|C10| ) yield equivalent results. The final answer contains only "point" 
mass terms, and is equivalent to combining the completely s-independent terms from the 
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derivatives of Q with only the "external" potential terms arising from any acceleration. The 
same approach holds for the terms G£, G\, and U^P^ 



2 • 

The vanishing of many potential contributions at order s° depends critically on the fact 
that the terms ultimately depend on the factor V fcZj V -1 , and that two of the indices are 
contracted, since only one index is allowed to be free. However, at 3PN order, this is no 
longer the case. A simple example is provided by a 3PN term proportional to v ■ J i> k U' l P$ k . 
Integrating over body #1, the combination of U£ xt with P^u) gives a contribution to the 
equation of motion 



(f) n <jk> 




—v{v k [ p*d 3 xm 2 y^-V Q V 1 (r~ 1 ) (\p 2 (f)5 jk + PT 2 ( 
m 1 Ji g q\ \6 

1 

p*f 2 PT 2 (f)d 3 x/m 3 mjm 2 v{v k V ijk (r~ 1 ) , (Cll) 



where the quantity in braces is dimensionless, scales as s° for fixed mi, but depends on the 
internal structure of body #1. Contributions like this appear everywhere at 3PN order; 
whether they survive in the final equation of motion, or what their ultimate interpretation 
is, will be the subject of future work. 



4. The quadrangle potential H 

The potential H = P(U' tj P 2 ) is an example of a more complicated "quadrangle" poten- 
tial whose fundamental ingredient depends on the field point and three source points. To 
see this, we write 

d 3 



H = ±- / J^U^(x')P^(x') 
An Jm x-x' v ' 1 v ' 



1 1 -(■■'■'{,■')— [ f X " U'\x")U>i{x") 



An Jm |x — x'| 4n Jm |x' — x" 

= E / P A V A V j A d 3 x A ( p B V l B d 3 x B f p* c V 3 c d 3 x c H(xA-BC), (C12) 

ABC Jb Jc 

where the function Ji of four field points is defined by 

H{AB] CD) = jk?U M \x A ~ x'| |x fl - x'| |x' |xc - x»\\x D - x"| ' (C13) 
with the properties 

V 2 A H(AB; CD) = -{l/r AB )Q{ACD) , 
V 2 C H(AB;CD) = -(l/r CD )G(ABC) , 
V A V B H(AB; CD) = (l/2r AB ) [Q(ACD) + Q(BCD)} - J(ABCD)/2 , (C14) 

where 

J(ABCD) = — / r . (C15) 

4tt Jm Ix-a — x' xb — x' xc — x' x D — x' 
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Unfortunately, we have been unable to find a closed-form solution for 7i similar to that for 
Q, nor a useful expansion in the case where some of the distances between points are small 
compared to the others. 

Instead, we make use of the first form of H given in Eq. (|C12j ). We integrate H' 1 over 



the density of body #1 and substitute U' k = U{ k + u£ and P{ k = P^ 11} + 2P 2 ( g 2 ) ) + P^ 22) . 
The result can be put into the form 

/ p'H'Wx = ^J m <Px'Ui{x>){Uf{x>) + Uf(x'))(P^ n) + 2P 2 fl + P§ 2) ) . (C16) 

We wish to verify that no contributions of order s° (other than normal point mass terms) 
arise in this integral. To see this, we split the integral over M. into an integral over body 
$T out to a radius ~ s%, a similar integral over body #2 to a radius ~ S2, and an integral 
over the rest of A4. In the latter integral, we may use solutions external to each spherical 
body: U = mi/j/i + m 2 /?/2 and P% from Eq. (pH). If carried out over all of M. using these 
functions, the integral would diverge at the locations of the two bodies Q . 



Consider now the integral over a region of volume s\ surrounding body #1. In the 
neighborhood of #1, the product UiUf k behaves as s J s J s fc /s 8 ; multiplying by the volume, 
we have a term of odd parity and 0(sj~ 2 ). The term Pjvii) ^ s even parity, so the combination 
integrates to zero. The term P^nl) mus t be expanded about #1 using Eq. (|U9|); the 



2(12) 

expansion begins at 0(s°) with a constant term and an odd parity term proportional to 
s l /s; then at O(s) with a term proportional to s (even parity) and one proportional to s l 
(odd parity); then at 0(s 2 ) with terms proportional to ss l and s l s J , and so on. All even 
parity terms integrate to zero when multiplied by s l s : >s k /s 8 . The odd-parity 0(s°) and 
0(s) terms lead to non-zero integrals of order sj~ 2 and sj" , which we discard. The odd- 
parity contribution of order s 2 , is accompanied by a coefficient proportional to V jfc '(l/r). 
Integrating over the sphere then results in a term of order s° but proportional to V 2 (l/r), 
which vanishes. Finally expanding the term P^l) about the location of body #1 gives only 
terms of order s m and parity (— l) m , hence only a contribution of order sj -1 survives in the 
integral. Applying similar considerations to the product UiU£ k , and then repeating the 
considerations for the integral over the neighborhood of body #2 leads to the conclusion 
that the contributions are of order s^ 2 or s^ 1 , or of positive powers, but that there are no 
structure-dependent contributions of order s° A . 

Consequently, we can carry out the integral in Eq. flC16|) over M. up to spheres sur- 



rounding each body, and then let the spheres shrink to zero, discarding terms that blow up 
as s^ 2 or s^ 1 , and keeping only finite terms. We are guaranteed that no structure-dependent 
terms of 0(s°) will appear. The final result is given in Eq. Q3.9|). 



APPENDIX D: EVALUATION OF MULTIPOLE MOMENTS FOR TWO-BODY 

SYSTEMS 

1. Quadrupole and higher moments 

Substituting expressions for from Eqs. ( |A5|) into the definitions of the multipole 
moments, Eqs. (|A4|) , converting from densities a, a % and <r y to the conserved baryon 
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density p* to the needed PN order using Eqs. (|2.21|) , and integrating, discarding any terms 
that depend explicitly on the radius 1Z of the boundary of Ai, we obtain 



A 1 AB r AB \2 4 J 

A 

J ij = e iab J2m A v b A x a /, 

A 

jijk... = e ^bJ2m A v b A x a f- , 



1 AR r AB 



-n A l B 5 ij + nf B 6 l)j + nf B 5 l)i - 28^ k 6 l)j + 2S ij 6 kl . (Dl) 



Note that, although X ijk and J ij appear in 2.5PN terms, they are purely functions of time, 
and thus cancel out of the relative equation of motion, so they are only needed to lowest 
order for use in 3.5PN terms. 

Converting to relative coordinates, using the 1PN correct transformation in the leading 



term of X %3 , we obtain 





= rar\x %3 (l + -(1 


— 3r])v 2 — 


J Tftj \ 7 

-(I-277)— J + -m 2 vr5 ij , 


j-ijk 


= -5mr]x ljk , 






j-ijkl 


= mrj(l — 3r])x lJ 






jij 


= —bmr\L % x 3 , 






j-ijk 


= 777,77(1 — 3rf)L l x 3 


k 

) 




M ijkl 


= 777,77(1 — 377) \v iJ 


1 777, ,-, N 

n 3 

3 r / 


)x kl 




_I m 2^ r (^ijfikl 


+ n kl 5 ij - 


n i{k 6 i)j _ n j{k 5 i)i + 25 i{k 5 l)j 



(D2) 



where L = x x v is the orbital angular momentum per unit mass. Note that the moments 
jijkim^ jijki anc j j^ijkim are no ^. neec i ec i explicitly for the 3.5PN equations of motion, since 

they are purely functions of time, and cancel out of the relative equation. 

Time derivatives of the moments may be calculated using the relative equations of motion 
in place of x % \ 1PN equations must be used in the leading term in X %3 \ while Newtonian 
equations are sufficient for the remaining terms. 



2. Dipole moment and the system center of mass 

For a two-body system, the dipole moment is given by X 1 = J M T 00 x i d 3 x. Substituting 
for r 00 including 1PN and 2.5PN terms from Eq. (5.9) of Paper I, and convering from a to 
p* to the corresponding order from Eq. ( 2.21D , we obtain, to zeroth, PN and 2.5PN order, 
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2 



(3) (3) ' 



-- E I <^ l kk +2 ^ I • (D3) 

Choosing the center of mass so that X % = to at least PN order, we see that the final 2.5PN 
term is in fact at least 3.5PN order. However, we must also check the time dependence of 
X 2 , to see if it remains zero to the order needed. From the definition of T 1 (Paper I, Eq. 
(4.6)), we have 

V = p l - / T 0j x i d 2 S j . (D4) 

JdM 

Using the definition of r - 7 (Paper I, Eqs. (4.4)) and the far- zone forms of the gravitational 
potentials (Paper I, Eqs. (5.12), with X % = to 1PN order), it can be shown that the surface 
integral is of 2.5PN order relative to P l . However, taking an additional time derivative and 
using T % i in the far zone gives 

V = P i = - 




(D5) 

where, to the order needed, X is simply the total baryon mass m of the system. Hence 
integrating with respect time and setting the initial conditions Xq = and Pq = 0, we have 

/ o (3) o (2) \ 
r = -m — V jj — e qik J qk . (D6) 

V 15 3 ) 

Note that, while this may seem like an anomalous 2.5PN effect on the system center of mass, 
it is really a gauge effect. Because the right hand side of Eq. ( P6|) is a total time derivative, 
it can be absorbed into a redefinition of spatial coordinates. Combining Eqs. ( p3| ) and 
flD6|), and defining x\ = (m-ilvn)x % + Q ', x\ = —im\jm)x % + Q ', we can then show that the 
transformation from x\ to relative coordinates is given by Eqs. ( |5.3| ), which were derived 
directly from the 1PN and 2.5PN equations of motion. 
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FIGURES 




FIG. 1. Transformation to the local comoving frame of one body. Global harmonic coordinates 
are t,x J ; local coordinates are t,x J . In local coordinates, the density of the body is assumed to be 
static and spherically symmetric. 
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